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1. Introduction 

Wind-driven sand transport is the most noticeable process shaping the morphology 
of arid regions on Earth, Mars and elsewhere. It is responsible for the spontaneous 
creation of a whole hierarchy of self-organized dynamic structures from ripples over 
isolated dunes to devastating fields of shifting sands. It also contributes considerably 
to dust proliferation, which is a major determinant of our global climate. There 
is thus an urgent need for mathematical models that can efficiently and accurately 
predict aeolian sand fluxes. But the task is very much complicated by the complex 
turbulent flow of the driving medium (e.g. streaming air) and the erratic nature of 
the grain hopping it excites [T). Yet, it is by now well understood that this hopping of 
grains accelerated by the wind has some characteristic structure [2]. Highly energetic 
"saltating" grains make a dominant contribution to the overall mass transport. When 
impacting on the sand bed, they dissipate some of their energy in a complex process 
called splash, ejecting a cloud of "reptating" grains [3111]. Grain-scale studies and 
theoretical considerations have indicated that a reduced description in terms of only 
two idealized populations (sometimes called saltons and reptons) should indeed be able 
to provide a faithful parametrization of the complex aeolian sand transport process 
and the ensuing structure formation [2j [5] . 

A drawback of the two-species description has been that it is still conceptually 
and computationally quite demanding. For reasons of simplicity and computational 
efficiency, many theoretical studies have therefore chosen to reduce the mathematical 
description even further, to mean-field type "single-trajectory" models This may 
be justified for certain purposes, e.g. for the mathematical modelling of aeolian sand 
dunes, which are orders of magnitude larger than the characteristic length scales 
involved in the saltation process, and thus not expected to be very sensitive to the 
details on the grain scale. Their formation and migration is thought to predominantly 
depend on large-scale features of the wind field and of the saltation flux, chiefly the 
symmetry breaking of the turbulent flow over the dune and the delayed reaction of the 
sand transport to the wind [7j. Moreover, the reptating grains, though they are many, 
are generally thought to contribute less to the overall sand transport, because they 
have short trajectories and quickly get trapped in the bed again. Therefore, it seems 
admissible to concentrate on the saltating particles. On this basis, numerically efficient 
models for one effective grain species that can largely be identified with the saltating 
grain fraction, have been constructed. The Sauermann model [8] is a popular and 
widely used example for such mean-field continuum models. One-species models have 
become a very successful means to gain analytical insight into EH EQ HI E3 Efl 
and to conduct efficient large-scale numerical simulations of [HI EH EO E3 ED. the 
complex structure formation processes caused by aeolian transport. The reduction to 
a single representative trajectory makes the one-species models analytically tractable 
and computationally efficient. But it is also responsible for some weaknesses both 
concerning the way how the two species are subsumed into one [5] and how they 
feed back onto the wind |I2j . These entail imperfections in the model predictions, 
most noticeable a systematic overestimation of the stationary flux at high wind 
speed (cf. figure [5j below). Therefore the one-species models have been criticized 
for their lack of numerical accuracy and internal consistency [5]. There are also a 
number of interesting phenomena that cannot be quantitatively modelled within a 
one-species model, because they specifically depend on one of the two species, or on 
their interaction, or because the two species are not merely conceptually but physically 
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different. 

As three important representatives of such phenomena, we name dust emission 
from the desert, ripples, and megaripples. Dust is created, exposed to the wind, and 
emitted by aeolian sand transport, and saltating and reptating particles play quite 
different roles in these processes [H E01 HI]- The reptating grain fraction and its 
sensitivity to the local slope of the sand bed is held responsible for the spontaneous 
evolution of aeolian ripples [H [4j E21 EJ EH EH EH ESI EH [28]. But the amount 
of reptating particles depends on the strength of the splash caused by the saltating 
grains. This interdependence of the reptating and saltating species becomes even 
more transparent if the two species correspond to visibly different types of grains. 
This is the case for megaripples, which have wavelengths in the meter range and form 
on strongly polydisperse sand beds in a process accompanied by a pronounced grain 
sorting [29 , 3Q [31] ■ The mechanism underlying their formation and evolution is still 
not entirely understood, but it seems that the highly energetic bombardment of small 
saltating grains drives the creep of the larger grains armouring the ripple crests. 

The demand for more faithful descriptions of aeolian sand transport has recently 
spurred notable efforts to eliminate the deficiencies of the one-species models [12] HU 
152] or to develop a numerically efficient two-species model They also motivated 
the development of the model described in detail below, as we felt the problem had 
not yet been cured at its root. In our opinion, many of the amendments proposed 
so far have targeted the symptoms rather than the disease by invoking some ad hoc 
assumptions. It was our aim to modify the sand transport equations from bottom-up, 
based on a careful analysis of the physics on the grain-scale and including the feedback 
of the two dissimilar grain populations on the turbulent flow. In this way, we derived 
a conceptually simple, analytically tractable, and numerically efficient two-species 
model with only two phenomenological fit parameters. They serve to parametrize the 
complicated splash process and take very reasonable values if the predicted flux law 
is fitted to measured data. Compared to the time when the first continuum sand 
transport models were formulated, we could build on a more detailed experimental 
and theoretical knowledge of the grain-scale physics [5] [331 13H [33] [SH] [37] [35] and rely 
on more comprehensive empirical information about the wind shear stress and sand 
transport to test our predictions [3j2 SO] HJ El] • 

The plan of the paper is as follows. We first summarize some of the pertinent 
basic notions of turbulent flows and aeolian sand transport and introduce the two- 
species formalism. Then we implement the two-species physics on the level of the 
sand flux in section 13.11 and on the level of the turbulent closure for the wind in 



section |3.2| In sections 3.3 and |3.4[ we address the most interesting observables, 
namely the speed and average number of hopping grains, and the saturated sand flux. 
We finally compare our results to other models and experiments in section [4] before 
concluding in section [5] 



2. Basic notions of aeolian sand transport and the two-species 
parametrization 

Our analysis of the two-phase flow of air and sand is similar, in spirit, to the Sauermann 
model [8 . The Sauermann model is a continuum or hydrodynamic model. Rather 
than with the positions and velocities of individual grains it deals with spatio-temporal 
fields of these quantities, viz. the local mass density g(x, t) and sand transport velocity 
v(x,t). It is of mean-field type since it reduces both fields to the mean quantities 
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Figure 1. The two-species picture of aeolian sand transport maps the 
ensemble of mobile grains onto two representative species. High-energy 
("saltating") grains accelerated by the wind (arrows in the wind-speed profile 
indicating the wind velocity) rebound upon impact and eject some low-energy 
("reptating") grains in a splash. Sand transport is quantified in terms of the 
transport velocities, i> sal and v Tep , and mean densities, p sal and p rep , of the two 
species. The momentum balance of the splash process (inset, cf. section \3.1\ is 
effectively encoded in the impact velocity u sal , the horizontal and vertical rebound 
velocities of the saltating grain, i)g al and v b Jq, respectively, and the number N and 
velocity u^ p = v™ p of the ejected reptating grains (neglecting a small horizontal 
velocity component, for simplicity). 



p(x, t), v(x, t) for a single representative trajectory, characterizing the conditions along 
the wind direction (the cc-direction), with the distribution of airborne grains in the 
vertical direction (the z-direction) integrated out. Note that p therefore is an area 
density, as opposed to the spatial volume density g. (The vertical grain distribution 
is a crucial element of the modified turbulent closure, though, which accounts for the 
feedback of the sand onto the wind.) In this contribution, we deal exclusively with 
the saturated sand flux, i.e. the sand flux along the wind direction over a flat sand 
bed and under stationary wind conditions. What we call the flux q — pv is actually a 
vertically integrated flux, or a sand transport rate [8]. The central aim is to derive a 
constitutive relation q(r), giving the stationary flux at a given wind strength r. 

In the following, the strong mean-field approximation of the one-species models 
is somewhat relaxed and reptating and saltating particles with densities p rcp , p sal and 
transport velocities u rcp , w sal are treated separately. The flux is split up accordingly 

q = q lcp + ? sal = p sal « sal + p Icp v rcp . (1) 

Introducing the dimensionless mass fraction (p € [0, 1] of saltating grains relative to 
the total amount of mobile grains, the (integrated) densities of the mobile grains can 
be written in the form 

p sal = ipp and p rcp = (l-cp)p. (2) 

Accordingly, the flux balance and the relative contribution of reptating and saltating 
grains to the flux read 

q = pv = p [ipv sal + (1 - tp)v rcp ] and 

q rep/ qSa .l = ^-1 _ iycp/ w sal _ ! - ,! 

The fluxes depend on the wind velocity field, which, in turn, depends on the 
grain- wind feedback. To make progress in this respect, the (steady) shear stress t in 
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the saltation layer is split up into two components carried by the airborne grains ( "g" ) 
and the air ( "a" ) itself, respectively, 

t = T g (z) + r a (z) = const. (4) 

This idea, introduced by Owen [43], which amounts to an effective two-fluid 
representation of the wind and the airborne sand, was used in many analytical models 
@1| El HUGS, and confirmed by numerical simulations H3 [Ml 133 [38] ■ In the two- 
species model, we moreover keep track of the two species of mobile grains in Q. 
Similar to the flux division introduced in ([!]) we further split up the grain-borne shear 
stress 

T g (z)=T^(z)+Tf(z). (5) 

Concerning the feedback of the grains on the wind, it turns out that the wind velocity 
profile u(z) is predominantly determined by the reptating grain fraction, and thus 
by the functional form of Tg' Gp (z), while the number and trajectories of the saltating 
grains hardly affect the wind profile, at all. The wind profile may thus be computed 



in a reduced one-species scheme, as explained in section 3.2 and Appendix C 

Under steady conditions, a common simplifying assumption is that the airborne 
stress near the ground (z — 0) can be identified with the threshold shear stress r t 
required to mobilize grains from the ground [43) . 

T a (0) PS T t . (6) 

If it were larger or smaller, an increasing or decreasing number of grains would be 
mobilized, respectively, resulting in an unsteady flow. This idea is indeed supported 
by empirical observations, which find that the number of ejected grains increases 
(almost linearly) with the impact speed, not their ejection velocity or jump height 
[4l)l |4"T1 1331 [31] . The feedback of the grains on the wind thus essentially fixes the air 
shear stress at the ground (and, if r m r t , also above) to the threshold value r t . To be 
precise, r t is called the impact threshold, to emphasize that it is easier to lift grains 
from the splash cloud than to lift completely immobile grains from the ground [T]. 
But, for the sake of simplicity, we do not bother to make this distinction, here, nor do 
we speculate about possible small deviations of r a (0) from r t under steady conditions 
[El [32]. Combining @ and we may express the shear stress contributed by the 
grains near the ground as 

T g (0)=TfP(0)+r! al (0) = T-T t . (7) 

Per definition, both stress contributions are proportional to the corresponding grain 
density. Following Sauermann et al [8] in assuming parabolic grain trajectories with 
initial horizontal and vertical velocity components v and v z0 , respectively, we find 

= gp/(2v zQ ) , r g (0) = gpiv - v )/{2v z0 ) , (8) 

for both species (i.e. species indices "rep" or "sal" may be attached to all variables). 
As a slight technical simplification, to which the following discussion is not sensitive, 
we choose to distinguish only the horizontal and vertical rebound velocity Vq &1 and 
v z f of the saltating grains (figure [l]) . Since the impact angle of saltating grains is 
known to be almost independent of the wind strength and to take typical mean values 
about 10° [48], we can identify the total impact speed with the horizontal speed i> sal 
of the impacting particle. Similarly, since the average ejection angle of the reptating 
particles is known to be about 80° [33J, we identify the total ejection speed of the 
reptating grains with its vertical component v z c ^. 
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In order to derive the sought-after constitutive equation, we have to understand 
the balance of the two species of mobile grains. If one considers the vertical fluxes 
sal and rcp rather than the horizontal fluxes <? sal and q lcp (the former follow from 
the latter via division with the jump length of the grains), one has 

^ c p = jV</> sal , (9) 

with the average number N of reptating grains ejected by an impacting saltating grain. 
Together with the first relation in Q, we then have 

ip- 1 =l + Nv T J/v^ (10) 

Using also the second relation in |8]) and remembering that we dismissed the horizontal 
component of the ejection velocity, we can thus write the overall density of mobile 
grains in the form 

_ 2g a j r U t — M„ t (Tl) 

P ~ g (p{v sal -v s al )/vll l + (l-tp)v™P/vl c p ' ( ' 

As often found in the literature, we have employed the notion of the shear velocity or 
friction velocity u*, defined by r = p a irU*, as a more intuitive measure of the shear 
stress r, here. The velocity ratios in the denominator of (11 1 are supposedly only very 
weakly dependent on the wind strength. The first one is recognized as the inverse of 
the effective restitution coefficient 

a = </(« 8al - ^ al ) (12) 

of the saltating grains [8]. Because of the complexity of the splash process, which 
makes first-principle quantitative estimates of a forbiddingly complex, we suggest to 
treat it as a phenomenological fit parameter. We do however provide a reasonable 
theoretical estimate of the ratio v rcp / w zo P ' based on the trajectory of a vertically 
ejected grain driven by the wind, in section 3.3 and Appendix D| Anticipating the 



result v Tep /v™Q = 0.7(cr^ a i r g) 1 / 3 /'u*t, we can rewrite the density of the transported 
sand in the more explicit form 

p _ 2gair u * ~ u *t (13) 

g ip/a + 0.7(1 - (^)(CT^air3) 1/3 /' lt *t 

Here u a i T denotes the kinematic viscosity of air, and, since the grain-air density ratio 
<t> 1 under terrestrial conditions, we do not bother to distinguish a and a — 1, here 
and in the following. Note that on top of the explicit wind-strength dependence of p 
via the numerator, there is an implicit, yet undetermined one via (p. 

In the following, we develop the sought-after "second generation" transport law 
q(u*), based on the grain scale physics. Some empirical input is employed to fix 
certain phenomenological coefficients summarized in table El Our result (see table [T]) 
accounts for essential elements of the aeolian sand transport process beyond the single- 
trajectory approximation, and turns out to be in remarkable agreement with wind 
tunnel measurements. 



3. Implementing the two-species framework 

The two-species framework is implemented in three steps. First, we deal with the mass 
balance between the species, then with their feedback onto the wind and their resulting 
transport velocities, before we finally arrive at an improved stationary transport law 
q{u*)- 
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3.1. Sand density: two- species mass balance 

We first estimate the partition of the grain population into a saltating and a reptating 
fraction from several evident assumptions for the splash process, based on fundamental 
physical principles. In our formalism, a saltating grain always rebounds upon impact 
and ejects N reptating grains, while a certain fraction of its kinetic energy and 
momentum is dissipated in the sand bed. We further assume that every reptating 
grain performs only a single hop after which it remains trapped in the sand bed. 
These mild approximations for the stationary sand transport reduce the number of 
free parameters of our model compared to more elaborate descriptions [S] . 

As argued e.g. by Kok and Renno [36], the relevant constraint limiting the 
number N of ejected grains per impact is given by the conservation of the grain-borne 
momentum rather than the energy, i.e. (cf. Fig. [T]) 

v ^ = v f + v l$ + Nv^ + v hcd (14) 

where the last term accounts for the bed losses. (Some omitted geometric factors 
will be subsumed in phcnomcnological numerical factors below, assuming a wind- 
strength independent scattering geometry). In a physically plausible linear-response 
approximation, all terms in the momentum balance are assumed to scale linearly 
in the velocity w sal of the impacting grain. For the first two terms on the right 



hand side of (14), this proportionality is already implicit in the definition (12 1 of 
the restitution coefficient a. There is an important subtlety concerning the last two 
terms, though. Namely, the ejected grains have to gain enough momentum to jump 
over neighbouring grains in the bed in order to be counted as mobile grains. But they 
also should not themselves eject other grains upon impact, otherwise we would have 
to count them among the saltating fraction. Hence, the impact velocity i> rcp and the 
ejection velocity i^g P of the reptating grains are tightly constrained and cannot, by 
definition, be strongly dependent on the impact velocity of the saltating grain or the 
wind strength. Within experimental errors, observations of the collision process of 
saltating grains by Rice et al [H] are indeed in reasonable agreement with a constant 
w^q P on the order of u*%. In fact, they found the horizontal component of the ejection 
velocity to be independent of w sal , and only a slight increase of the vertical component 
with w sal . A compelling confirmation of this observation was recently obtained by 
better controlled laboratory experiments, in which PVC beads were injected at an 
impact angle of 10° onto a quiescent bed of such beads [33J 131] ■ There are two 
important consequences of the constraints on v T z °^. Firstly, N should grow linearly 
with w sal . And secondly, since the momentum of the impacting grain is only partially 
transferred to the ejected grains, a critical impact velocity w sal = t^ al 3> vJcm nas *° 
be overcome to mobilize any grains at all. It can be interpreted as a constant offset in 
i>bcd, which does not scale with w sal . According to the collision experiments [331 134) . 
one can take i^ 3,1 « 40^/gd » 10u*f (The final expression for w^ al in terms of the 
aerodynamic entrainment threshold u*t is based on Bagnold's estimate u* t ~ Q.Qlagd 
PQ, with the PVC-air density ratio a = 1444.) A concise summary of the empirical 
input entering our model is provided in table [ET] 



Summarizing the above discussion, we can rewrite the momentum balance ( 14 ) 

as 

Nv^cxv^-vf, (15) 

where v sal > uj! al is tacitly assumed to hold throughout our discussion, and the omitted 
factor of proportionality should be insensitive to the wind strength. Observations of 
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Figure 2. The wind speed u at height z. Theory, equation \21\ (curves), 
compared to measured data (symbols) for grains of diameter d = 242 fim (solid 
curves, circles) and d = 320 (dashed curves, boxes) and various shear velocities 
(u* = 0.27, 0.39, 0.56, 0.69ms" 1 and u, = 0.27, 0.47, 0.74, 0.87 ms" 1 from 
bottom to top for d = 242 /im I41| and d = 320 fira 1401 . respectively); as global 
fit parameters we used the roughness length (20 = d/20), the threshold shear 
velocity («*t = 0.19 ms — 1 , 0.20 ms -1 for the finer/coarser sand), and the mean 
saltation height z m , which is compared to \23\ in the inset. 



saltating grains [15] , model collision experiments [15J EH |M] , particle dynamics 
simulations and reduced numerical models, such as the binary collision scheme 

proposed by Valance and Crassous |35) . all support this afhne relation of the number 
of ejected grains, which is the cornerstone of our two-species mass balance. To fix the 



omitted factor in ( 15 ) and make contact with ( 10 ) we write 



Nv I ;» = (r,/a)(l-v^ /« sal )< . (16) 

The constant 77, which serves as the second free fit parameter of our model, determines 
together with u^ al how the momentum lost by a rebounding grain is distributed 
between the bed and the ejected particles. (Details concerning the phenomcnological 
values of a and 77 obtained from fitting the complete model to experimental flux data 



are given below.) Comparison with (10) now yields an explicit expression of the mass 
fraction 

p" 1 = 1 + (1 - vf/v^/a (17) 

in terms of v sal , alone. Thereby, the problem of specifying the stationary mass balance 
between saltating and reptating grains has completely been reduced to the task of 
finding the stationary transport velocity v sal of the saltating particles as a function of 
the wind strength. 



3.2. Wind velocity field: two-species turbulent closure 

Before we can work out the transport velocities of the two species of grains, we need 
to know the height-dependent wind speed u(z), and how it is affected by the presence 
of the airborne grains. Unfortunately, this leads us back to the question of the grain 
densities that we just delegated to the calculation of the grain velocities. In the past, 
an elaborate self-consistent calculation has been avoided by anticipating the resulting 
form of the wind profile on the basis of empirical observations. In the Sauermann 
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model, for instance, an exponential vertical decay of the grain-borne shear stress is 
imposed, 

r g (z) = T g (0) exp(-z/z m ) , (18) 

in good agreement with grain-scale simulations (6j |45j [36l [37]. Since the airborne 
stress r a (z) follows from this via Q, it reduces the task to the problem of finding 
the mean saltation height z m . In [Appendix C| we present a refined version of this 
approach, adapted to the two-species approach. It accounts for the contributions of 
the different species to the grain-borne shear stress, separately, and also for their 
considerably different trajectories. However, as we detail in the appendix, the whole 



exercise yields essentially identical results as (18 1, which can in fact be interpreted 



as the pre-averaged two-species expression. This lends further support to the general 



form of (18 1, recommending it as a suitable basis for our further discussion of the 
two-species model. 

From (18), we get, via Q, the wind speed profile u(z) by integrating Prandtl's 



turbulence closure 

T a (z) = Q^ 2 Z 2 [d z u{z)] 2 ■ (19) 

To proceed analytically, we use a modified secant approximation similar to what was 



proposed by S0rensen [11] . As shown in Appendix B , the closure equation then 
becomes 



KZ 



1 - (1 - u* t /u*)e~ z/Zm . (20) 



Upon integration from the roughness length z <c z m , defined by u(z Q ) = 0, this yields 
the explicit result [TTJ [T2] 

u{z) = ^ m(z/zo ) _ u * - "** [ El ( Zo / Zm ) - El (z/z m )} , (21) 

K K 

with the integral exponential Ei(z) = J°°dx e~ x /x. For small wind speeds — > it*t 
(u* > U*t), the usual logarithmic velocity profile is recovered. Inside the saltation 
layer, a universal asymptotic wind velocity field 

u(z) ~ ^ ln(z/z ) (z < z « z m ) (22) 

K 

emerges, which is independent of the wind strength outside the saltation layer. This 
rationalizes the convergence of the wind profiles for different wind strengths, found in 
wind tunnel measurements, without the ad hoc assumption [12] of a debatable focal 
point _50i, which resulted in the prediction of an unphysical negative flux q(u*) at 
high wind speeds |12j . A direct comparison of our above prediction for u(z) with wind 
tunnel data from references [40] and [41] can be found in figure [2] From these data 
(available for two grain sizes d) we determine the mean saltation height z m . Past 
attempts to relate z m to d invoked a new undetermined length scale depending solely 
on atmospheric properties [ST]. However, the results of our fit rather support the 
simpler (and more natural) linear relation 

z m « 25d , (23) 

(right panel of figure [2]). This is in line with recent analytical and numerical work 
[32j , albeit possibly with a somewhat different numerical factor dependent on specific 
model assumptions. 
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Figure 3. The transport velocities of the saltating (dashed) and 
reptating (dotted) sand fraction versus the rescaled shear velocity 

U = u,/u,t, as predicted by l |24[ ) and l |27| , respectively. The solid line represents 
the average transport velocity in a reduced effective single-trajectory description. 



3.3. Transport kinetics: two-species transport velocities 

With the wind speed profile u(z) at hand, we can now determine the characteristic 
stationary transport velocity w sal of the saltating sand fraction from a standard force- 
balance argument. We evaluate the wind speed at a characteristic height z sal , at 
which the fluid drag on the grains is then balanced with the effective bed friction, in 
the usual way [12] . Namely, the drag force on a volume element of the saltation 
cloud (the force per mass acting on a single saltating grain times the density p sal of 
the saltation cloud) is balanced with T| al (0). This yields the transport velocity of the 
saltating sand fraction 

y sal = u ( z sal ) - v ^ , (24) 



with the estimate (52l [12 



t£ « ^ Jl9la r (25) 



for the terminal steady-state relative velocity u — w sal . Formally w^ 1 is equal to the 
settling velocity asymptotically reached by a grain freely falling in air (of kinematic 
viscosity f a ; r and with grain-air density ratio er), with a rescaled gravitational 
acceleration g M- gj (2a). 

We fix the -independent parameter z sal by the plausible assumption that 
the splash at the (impact) threshold wind speed it*t dies out, corresponding to 
w sal («* = u*t) = v^ 3,1 . The underlying physical picture is that it is the splash that 
keeps the saltation process going, even below the aerodynamic entrainment threshold, 
because the reptating particles compensate for rebound failures [5] . Accordingly, using 



(24 1 with the logarithmic wind field, equation (22 1, we find 

sal i ^.sal * 



z sal = z exp( k — °° ) . (26) 



it-* 



In |Appendix D[ we use a similar argument to derive the mean reptation velocity, which 
turns out to be almost independent of the grain size d. It can thus be approximated 
by the constant 

W r °P w 0.7(a^ airff ) 1 / 3 , (27) 
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for most practical purposes. This relation corresponds to a -independent reptation 
length 2i> rep uJg P / 1 g = 0.14[cr 5 zA r <i 3 /g] 1 ^ 6 in the centimeter range, while the saltation 
length 2v sal v s z f/g oc v%/g is quadratic in the shear velocity (cf. figure [3]) and on 
the order of decimeters [TJ [531 HI] ■ 

In figure [3J the transport velocities u sal and t> rcp for the two species and the mean 
velocity v = yw sal + (1 — (p)v rcp are plotted over the shear velocity u*. An important 
result is that the species-averaged transport velocity v is nearly constant over a broad 
range of wind strengths. This is consistent with the fundamental assumption on which 
the model is based, namely that the number of mobilized grains is sensitive to the wind 
strength, but their overall transport kinetics is not. 



3. 4- Stationary sand flux: two-species transport law 

The stationary flux law is now obtained by collecting the above results for the saturated 



density p from (111, the saltation fraction <p from <jX0J) , the transport velocities (24) 
with the wind velocity (21), and inserting them into Q3J) . We introduce reduced 
variables, measuring the shear velocity w* in terms of the impact threshold u*t and 
the saturated flux q in units of Qa,hU*/g, 

U=u*/u* t Q = qg/ie^ul) . (28) 

The result for our two-species stationary flux relation then takes the form 

Q=(l-U->) a ° U+ ^:^ U ~\ (29) 

with the coefficients a 0l /3 , 70, o■^ and b depending on the two free parameters a and 
77, as summarized in |Appcndix E| It goes without saying that this result pertains to 
U > 1 and that Q(U < 1) = 0. Except for b, which becomes negative for small ?y, and 
in particular for 77 = 0, all coefficients are positive. But since 70 < ao + A) and b < a, 
the flux always remains positive. 



4. Discussion 



Before we determine the free parameters a and 77 by comparing ( 29 1 to experimental 
data, we first discuss the effect of our two-species parametrization on the two related 
transport modes — saltation and reptation — and present two different single trajectory 
reductions of the two-species model. Subsequently, we compare these reduced schemes 
to single-species transport laws that have previously been proposed in the literature. 



The two- species flux balance 



Our derivation of the overall sand flux Q, given in (29), also yields the partial fluxes 
gsai = q sal g / (g a i T ul) and Q rcp = q Icp g/{ea.irul) of the saltating and reptating species, 
respectively. Their dependence on the wind velocity is illustrated in figure |4j It 
is frequently argued in the literature that the saltating grains dominate the sand 
transport, for they perform large jumps, while the reptating sand fraction might be 
negligible in terms of mass transport. Our model basically confirms this view for 
moderate winds, not too far above the impact threshold. However, it predicts that 
both species contribute almost equally to the total flux, under strong winds. The shift 
in the mass balance ip towards reptation (figure |4j right plot) reflects the dependence 
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Figure 4. The rescaled fluxes Q sal = <? sal g/(0air"2) and Q rcp = 
<? rep 9/(?air u *) of the saltating and reptating grains, respectively, versus the 
rescaled s hear velocity U = u r ./u*t (parameters: a = 0.65 , r) = 10, as determined 
in section [4.3| grain size d = 250 (im.). Left: direct comparison of Q sal (solid 
line) and Q rc P (dashed line). Right: The flux ratio Q rc P/Q sal = g re P/g sal (solid 
line) compared to the mass fraction 1 — ip of reptating grains (dashed line). 



of the number of ejected grains on the wind strength, which follows from the splash- 
impact relation ( 15 ) if one inserts the transport velocities of the individual species, 
as given in figure [3j The variable contribution of the two species to the overall sand 
flux is a result that cannot be obtained from the conventional single-species models, 
but should be of interest for some of the more advanced applications mentioned in the 
introduction. 



4-2. One-species limits 

It is interesting to see how the conventional one-species descriptions emerge from the 
two-species model upon contraction to a single effective grain species. This contraction 
is clearly not unique but may be performed in different ways. On the basis of the 
two-species flux balance, discussed in the preceding section, two natural contractions 
suggest themselves: one for moderate wind speeds where the flux is dominated by 
saltation, and one for strong winds, where the ratio of the contributions from reptating 



Table 1. Scaling functions f(U) for stationary aeolian sand transport relations 
of the form Q(U) = (l — U~ 2 ) f(U). Note that the coefficients of the various 
models cannot be identified, even if represented by the same symbol. Explicit 
analytical expressions for the coefficients ao, Aj, 70i o,, and b of the two-species 
model can be found in | Appendix E| 

„, . . agU + Pa - loU' 1 

this work 

aU - b 

This work, one-species limit a — bll" 1 (b < 0) 

Sauermann et al 2001 [8] as/l + a JJ- 2 + bU' 1 

S0rensen 2004 [IT] a + bU' 1 + cU~ 2 

Duran and Herrmann 2006 [12] a- bU' 1 (a, b < 0) 

Pahtz et al 2011 [32] -a + bU' 1 + cU~ 2 
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Figure 5. The sta tionary sand-flux laws. Predictions from th e tw o-spe cies 
model, equation ^29\ (solid line) with its natural one-species limits ( | 30 [ > and | |31[ ) 
(lower and upper dotted line), are compared to transport laws from conventional 
one-species models: Sauermann et al [8] (upper dashed), Duran and Herrmann 
1 1 21 (lower dashed), and Pahtz et al 1321 (middle dashed). All curves have been 
obtained for a grain size of d = 250 fim. The values of the free parameters a and r) 
in ^29\ were determined by comparison to experiments as described in section |4,3| 
The predictions of the single-trajectory limits of the two-species model, \30\ and 
(| 3 1 p , were adjusted by the parameter rescaling z sal — > 0.18z aal and rj — t 0.138»7. 



and saltating grains was found to be roughly constant. 

The first single trajectory reduction, where the reptating species is dismissed, is 
obtained by setting tp = 1 or rj = 0, which corresponds to the flux relation 

Q = (1 - U- 2 ) (a- blf- 1 ) , (30) 



with a negative coefficient b < (Appendix E). Since this equation accounts only 



for high-energy saltating grains, the resulting absolute flux is too large compared to 
the original two-species description. This can be corrected by reducing the effective 
trajectory height or, in our formalism, by rescaling the effective height z sal at which 
the fluid drag is balanced with the effective bed friction. The lower dotted line in 
figure [5] was obtained in this way. 

The second single trajectory reduction is motivated by the weak dependence of 
the mass fraction ip on it* for strong winds it* 3> which suggests to replace tp by 
its -independent limit 1/(1 + Tj/a). This yields 

Q={l-U- 2 ) (<S + A) t/" 1 ) , (31) 

with 

"o , 3 A) . V v T z l P 
a = — and /3 = — + — — — - , (32) 
a 2a 1 + 77 u* t 

and the abbreviation 77 = r]v rcp / v 1 ^ . If one interprets the (now constant) species mass 
balance as a free phenomenological parameter, it may be adjusted by fine-tuning 77, so 
that the absolute value of the sand flux predicted by this formula agrees better with 
experimental observations. This is how we obtained the upper dotted line in figure [5] 



Both (30) and (31) have the functional form of the transport law proposed by 
Duran and Herrmann [12 based on the focal point assumption. Note that the effective 
height z sal is small compared to the saltation height, where the wind speed obtained 
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Figure 6. The rescaled saturated flux Q = gg/(g a ir u 2) f° r various grain 
sizes (d = 125 fim, 170 fim, 242 fim, 320 fim, 242 fim from bottom to top): data 
from wind tunnel measurements using particle tracking methods 42 (stars) and 
sand traps 39 ( othe r symbols) are compared to the prediction by the two-species 
model, equation j29| (solid curves). The restitution coefficient a served as a global 
but grain-size dependent fit parameter (see figure n\ t while rf = 10 and r\ = 4.2 
was fixed for the sand trap and the particle tracking data, respectively. (The open 
symbols represent several rescaled data sets, obtained for a variety of bed slopes.) 



from the focal point assumption is very weak. As a consequence, the mean drag force 
is reduced, which results in a negative parameter a < in ( |30[) and a negative sand 
transport rate Q < for large shear velocities U > \a/b\. Table [lj and figure[5| provide 
a summary of various stationary aeolian sand transport laws that have been proposed 
in the literature, in comparison with the prediction of our two-species model and its 
above single-trajectory contractions. 



4-3. Comparison to experiments 



We now fit our two-species flux relation Q(U) from (29) to empirical data from 
various wind tunnel measurements [391 1421 . using r\ and a as free fit parameters. As 
demonstrated in figure[6j we obtain excellent agreement with the experimental data for 
a grain size independent splash efficiency r\ and an effective coefficient of restitution a 
that increases with increasing grain size d. The transport rates for four different grain 
sizes collected in reference |39j were obtained using sand traps; the data in reference 
|42j was obtained by particle tracking techniques for a single grain size only. The 
different methods yield different absolute values for the sand flux, which might be due 
to systematic differences in the detection efficiency. In the model, these differences 
are reflected in (apparently) different efficiencies with which saltating grains generate 
reptating grains, i.e. different values of the parameter rj. We find rj = 4.2 for the data 
obtained by particle tracking and rj — 10 for the sand trap data, independently of 
the grain size. But both data sets were consistent with the same formula for the bed 
restitution coefficient 

a = 1 — do/d , (33) 

with d > do — 88 /im (left panel of figure [7]). That the rebound becomes increasingly 
elastic for large grains seems physically reasonable. Also that a vanishes for a critical 
grain diameter d = do is in accord with the observation that saltation only occurs for 
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Figure 7. The restitution coefficient a (left) and the threshold shear 
velocity u»t (fight) versus the grain size d. Parameter values deduced from 
model fits to flux data in figurejfi] The solid lines represent the scaling (d/dg)^ 1 
according to d 33 b for 1 — a and the classical estimate u*t = O.ly/crgd, respectively. 



sand grains larger than about 70 /jm |55j . The collision of smaller grains is increasingly 
influenced by hydrodynamic interactions (and besides by cohesive and electrostatic 
forces). In other words, do is a characteristic grain size marking the transition from a 
phenomenology typical of sand to one typical of dust. 

The observed value of do can be rationalized by a hydrodynamic order of 
magnitude estimate. The idea is that the crossover from sand-like behaviour to dust- 
like behaviour is concomitant with the crossover from turbulent to laminar flow on 
the grain scale. Sand grains collide with each other and with the bed, and they are 
lifted from the bed by turbulent fluctuations in the wind velocity. In contrast, dust 
particles are essentially hidden from collisions and from turbulent fluctuations trying 
to lift them from the bed by the laminar boundary layer that coats any solid (no slip) 
boundary. The characteristic length scale of this laminar coating is related to the so- 
called Kolmogorov dissipation scale (iA^e) 1 / 4 [55], where e is the dissipation per unit 
mass of the driving medium (in our case air at normal conditions) . For the logarithmic 



wind profile near the ground, equation (22), the scale dependent dissipation takes the 
form e(z) = u^ t / ' nz [57]. A self-consistency argument requiring that the Kolmogorov 
dissipation scale at height z above the ground equals z (if no externally imposed 
roughness scale is available), yields z ~ [^ ir /£(d )] 1 ^ 4 - Extrapolating the logarithmic 
profile with u* t = 0.1 \Jogd down to the scale z, we obtain z = 4.6(i/f ir /crgf) 1 / 3 . We find 
this estimate to be in excellent agreement with the data in the left panel of figure [7] 
which suggests 

d o «4(i&/<70) 1 / 3 «lO a /an. (34) 

Moreover, while we also fine-tune the threshold it*t for each grain size by hand 
when fitting the experimental data, the values used turn out to be in very good 
agreement with the expectation u*t ~ 0.1y/agd, as demonstrated by the inset of 
figure [6] Altogether, the stationary sand transport rate observed in experiments is 



thus convincingly reproduced by the two-species result, equation (29 1, over a wide 
parameter range, with very plausible and physically meaningful values for the model 
parameters. 
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5. Conclusions 

We have developed a "second generation" continuum description of aeolian sand 
transport, relaxing the strong mean-field approximation inherent in the classical single- 
trajectory models such as reference jS], and introducing a two-species framework, 
as advocated in reference [5]. We started from the physics on the grain scale and 
corroborated our explicit analytical expressions by a comprehensive comparison with 
empirical data for the splash, the wind velocity field, and the sand flux, covering a 
broad range of ambient conditions and grain sizes. While the usefulness of our general 
approach of contracting the complicated splash and saltation process into a drastically 
reduced two-species picture was strongly suggested by empirical and theoretical work, 
it necessarily involved some free phenomenological coefficients. For our two-species 
model, these are, apart from a few minor numerical coefficients listed in table |El} the 
effective bed restitution coefficient a(d) for the saltating grains of diameter d, and the 
parameter rj that fixes the number of reptating grains ejected by the average impacting 
saltating grain. Both parameters were found to take consistent and plausible values 
if used as free parameters when fitting empirical flux data obtained in wind tunnel 



experiments. The predicted two-species stationary flux law (29 1 was found to be in 
excellent agreement with comprehensive data from different sources. It should provide 
an excellent analytical starting point for a variety of advanced applications calling for a 
more faithful description of the saltation process as so far available — from wind-driven 
structure formation in the desert to saltation-driven dust production and emission. 
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Appendix B. Solving the Prandtl turbulent closure 

With Q and the assumption of an exponential decay of the grain-borne shear stress 



with height, equation (181, the modified Prandtl turbulence closure (19) reads 



d z u= — ^ e -*/*». (B.l) 

KZ V T 

To make analytical progress, we follow S0ren sen |llj in approximating the square root 



by a secant. The right hand side of (B.l) has the functional form ^/l — ex, with 



x = exp(— z/z m ) and e = r g (0)/r, which we approximate by a secant of the form 
ax + b(l — x). Matching the points {0, 1} and {1, y/l — e} yields 6=1 and a = y/l — e, 
hence 

y/l - ex « 1 - (1 - y/l -e)x. (B.2) 
Under stationary conditions, using (M), we have y/l — e = y/(r — T g (0))/r = u*t/u* 



which leads to ( 20 1 . This avoids an artifact of S0rensen's original approximation [TT] , 
y/l — ex s» 1 — ex, which implies ^/r a (0)/r = 1 — r g (0)/r for the shear stress at the 
ground, as already criticized by Duran and Herrmann [12 . 



A two-species continuum model for aeolian sand transport 



17 



Appendix C. The two-species approach for the wind profile 



In ([5]) of section [2j the grain-borne shear stress was split up into the contributions 
from saltating and reptating particles, respectively. The ratio of the grain-borne shear 
stresses at the ground, 

T g rC P(0) V rQ P l-<p__l-tp 



r sal 



(0) 



rep 



'■P 



9 



(C.l) 



immediately follows from i8l and (12 1, and from our simplifying assumption that the 



reptating grains are ejected vertically. Assuming the exponential decay of the grain- 
borne shear stress r g with height to hold for both components, the Prandtl turbulent 
closure reads 

- Z / Z „T _ T- sal (o)e~^»' 



PairK Z (dgU)" 



r rcp 



(OK 



(C.2) 



g V"/~ -g 

We exploit the strong scale separation z^ l /z™ p ~ 10 2 between the characteristic 
jump heights of saltating and reptating grains [TJ[5H], on which the two-species model 
is based. (The precise value turns out to be irrelevant to our discussion.) It allows to 
solve the closure for two separate height-ranges: (i) z < z™ p -C z%~ associated with 
reptation, where we may set exp(—z/s%j~) — > 1, and (ii) z ^> z™ p , associated with 
saltation, where we may set exp(— z/z™ p ) — > 0. Applying the secant approximation 
for the square root as described in | Appendix B[ we can perform the integrations within 
both ranges and match the asymptotic solutions at z 
eliminate r| al (0) and r rcp (0), this yields 



y rcp 



Using and ( |C.ip to 



u{z) 



5w* 



, z Su* 
In 



zo 



In 



K 

Z 



rep 
Zm 



'Ei(z /z 



Mz/z^)} 



u»(l-5) 



E 1 (z2 P /z^)-E 1 (z/z^)] 



z < z™ p 



z > zt 



(C.3) 



where we abbreviated 



s = Wi- 



ll - T t /r)<p 



tp + (1 - ip) a 



Using this in the force balance, (24), we gain the implicit equation 

sal „,/~sal „,sal\ 



u(z s 



) 



,,sal 



(C.4) 



(C.5) 



where the fraction tp of saltating grains depends itself on the velocity v of the 



saltating grains via ( 10 ). To solve it, we note that the fraction tp of saltating particles 
is a decreasing function of the shear stress r, bracketed by 1 and (figure [4j right 
panel). Hence, the product (1 — Tt/r)tp is a small number for all r > r t . Setting it to 



zero in (C.3), we get 

, v u* f z 
u(z) ~ — In — 

K \Z 



u*-u* t f [E!(z /C P ) 
« X I IMzo/z^) 



Ei(^/C)] 
Ei(l)], 



7 < 7 r<=p 

z > z™ p 



(C.6) 



Since the exponential integral Ei(z/z m ) vanishes rapidly for increasing z > z m , we 
recover the wind speed profile (21 1 successfully employed in the one-species models, 
with the characteristic decay height given by z m = z™ p . Inserting (C.6) into (C.5) 
yields an affine increase of v sal with u*, in good accord with the exact numerical 



solution of ( C.5[ ) (figure CI right panel) 
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Figure CI. Comparison of the wind velocity profile | |C,3| obtained from the self- 
consistent numerical solution of | |C,5| l (solid lines) with the pre-averaged profile 
| |21| (dashed lines) for grains of diameter d = 242 ^m and for various shear 
velocities (U = 1.5, 2.5, 3.5 and 4.5 from bottom to top). As explained in the 
main text, both approaches coincide for z m = zJJ p . Addi tiona lly, we show the 
results from numerical solutions of the turb ulence closure ( |C.2| (open symbols) 
to support the approximate expression (C.3l. Inset: The saltation velocity tj sal , 
rescaled by the minimum saltation velocity t^ al needed to eject any grains, over 
the shear velocity: self-consistent numerical solution of | |C.5| (symbols) and its 
analytical prediction in the limit (1— T%/r)(p — > 0, for which u(z) is given by JC.I 
(line) . 



For the numerical solution, one has to deal with the two free parameters z™ p and 
a. While the former is directly related to the wind profile, the latter is a fit parameter 
of the model, determined from a comparison of the predicted sand transport rate 
with empirical data (figure [6]) . To make progress, we vary z™ p and take the value 
of a from section |3.4[ where the sand flux is estimated by means of the pre-averaged 
approach for the wind profile. For a grain diameter d — 242 /an, we obtained a ~ 0.63 
(i.e. a = 1.4), which is consistent with collision experiments, as argued in reference 



12] . From the numerical solution of (C.5), we find a good agreement between the 



self-consistently gained u(z) and wind tunnel measurements [41] for z™ p = 25d, which 
is exactly the same value as obtained in section [3 . 2| within the pre-averaged approach. 



This supports our observation that (C.3) can be approximated by taking the limit 



(p — > in ( G.6 ), which yields the pre-averaged wind profile for z m = z™ p . This result 



is almost independent of the ratio z^ /z^ p , for ( C.6 1 is independent of z^ . Thereby, 



we formally confirm the intuitive expectation that the feedback of the grains on the 
wind profile is predominantly due to the many reptating particles, and hardly affected 
by the saltating particles. 

Appendix D. Reptation velocity 

The (horizontal) velocity of a reptating grain obeys the drag relation 

d tv r0p = t4p^ \u - v Icp \ (u - z/°P) , (D.l) 

similar to the one for saltating grains [H1[T2]. We assume that the ejection is essentially 
vertical with the initial velocity i^q P on the order of u„ t . A more accurate discussion 
would not substantially change our findings, as confirmed by the numerical solution of 
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Figure Dl. The Reptation velocity v vcp against the grain diameter d. The 
dots correspond to the numerical solution of (JdTTJ averaged over one trajectory, 
while the solid and dashed lines represent the approximate solution jD.3[ l and 
the representative value v Icp = 0.7(<ri' a i r <?) 1 / 3 0.5 ms -1 proposed in ^27} , 
respectively. The latter can be understood as an average over the relevant grain 
sizes (white background). 



(D.l ). Note that the nonlinearity of the friction law entails a time-dependent "terminal 
settling velocity" v™ p , dependent on the actual relative grain velocity u — v rcp |52j . 
However, for our purpose, and in view of the low reptation trajectories, we can safely 



approximate the reptation velocity from (D.l I by inserting the wind speed u(z rep ) at 
a given reptation height and the steady-state terminal velocity 



„ re P _ 



0.95 + 20^ air /v / cr3d 3 



(D.2) 



similar to (25). Neglecting moreover vertical drag forces, the maximum height of the 
reptation trajectory is {v™^) 2 / (2g) « Wd < z m . Consequently, we may insert the 
ground-level wind field ( 22 ) and obtain 

v r Cp ra u *t ln(z rc P/z ) 

~ k \l + k (vZ p f/ (2<o P ^ t ln(z-P/z ))l ' 



(D.3) 



where we approximated the time of flight by that for a parabolic trajectory, 2i^q P / 'g, 
as usual. Inserting the empirical observation v T z op ~ u* t a s well as v% t = O.Olagd and 

as a 



z = d/20, cf. table El results in an estimate for the reptation velocity w rcp 
function of the grain diameter d, illustrated by a solid line in figure |D1| 

The yet undetermined value of the effective reptation height z rep is expected to 
be on the order of the maximum height (v r ^) 2 / (2g) of the trajectory. Indeed, if we 
compare the approximate result given by ( |D.3[ ) with the numerical solution v lcp (t) of 
(D.l ) averaged over the whole trajectory, 

r2v lc f 

9 



v rcp 



1v icp 



*/9 



dtv iep (t), 



(D.4) 



we obtain a good match for large grain diameters d > 800 /im for z rcp = 0.14uJ t / (2g). 
Note, however, that the numerical solution yields an almost (i-independent reptation 
velocity v Tep ~ 0.5 ms" 1 for the relevant grain sizes d w 100 /im ... 1 mm, as illustrated 



in figure Dl To find an estimate for i> re P, which still captures the dependence on other 



parameters appearing in (D.3), we evaluate this equation for a representative grain 
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diameter within the relevant range. A closer look at (D.3| reveals that the reptation 
velocity at the inflection point d ~ [ag)~^^ s=s 150 jj,m can be approximated 

by i> rcp ~ (v 'air o -g) 1 ' ^ = 0.68 ms _1 , which provides the wanted parameter dependence. 
To better match the absolute values found from the numerical solution, we insert a 



factor 0.7 by hand, thus arriving at the analytical estimate given in (27). 



Appendix E. The coefficients of the transport law 

Here we give the explicit expressions for the coefficients occurring in the saturated 



sand flux ( 29 ) 



and 



a = 2a(l + fj) 

„,sal 



Zq exp k 



„.sal i .sal 
U*t 



b = a - 2c 



n\2 



(1 + 7?) 

2f)a 



1 + fj V u*t 



To 



~ rep 
VV z0 



~ rep 



1 + JJ 

- 2{a - b) 



with the abbreviations 

£ i z ) = Ei(zo/zm) ^ Ex(z/z m ) 

and 



1 + fj 



(E.l) 
(E.2) 

(E.3) 
(E.4) 

(E.5) 
(E.6) 

(E.7) 

(E.8) 



Table El. The Parameters occurring in equations | |E.1| to |E.5| which are 
either numerical constants or dependent on the grain diameter d, the gravitational 
acceleration g = 9.81 ms — , the sand-air density ratio cr = 2163, or the kinematic 
viscosity of air, i/ a ; r = 1.5 X 10 -5 m 2 s — 1 . 



parameter value or formula 
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All parameters occurring in these equations are listed in tabic [El| where, as explained 
in the main text, the reptation velocity v rcp is independent of d and it*, the critical 
impact velocity u® al and the vertical ejection speed v Y z e ^ are obtained from collision 
experiments, the roughness length zq and the height z m (the characteristic decay height 
of the air-borne sand density) , are estimated by fitting experimentally observed wind 
velocity profiles above the saltation layer, and i], u*t, and a are determined form 
fitting the sand transport rate to wind tunnel measurements. 
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